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Abstract 

We study, using dimer and Monte Carlo approaches, the critical properties and 
finite size effects of the Ising model on honeycomb lattices folded on the tetrahedron. 
We show that the main critical exponents are not affected by the presence of conical 
singularities. The finite size scaling of the position of the maxima of the specific 
heat does not match, however, with the scaling of the correlation length, and the 
thermodynamic limit is attained faster on the spherical surface than in corresponding 
lattices on the torus. 
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1 Introduction 



The Ising model conveys, in its simplicity, a richness of physical information which makes it relevant 
as a model for critical phenomena in different instances (ferromagnetic materials, lattice gas, binary 
alloys). The model is also paradigmatic of a common situation in statistical physics since, being one of 
the simplest models, yet it only allows to compute analytically the thermodynamic limit for particular 
classes of lattices (in one or two dimensions) . In two dimensions, the Ising model in a square lattice has 
been solved in the continuum limit with cylindrical and toroidal boundary conditions 0, ^ • It has been 
also found analytic solution for the two-dimensional model on a triangular or honeycomb lattice . In 
general, however, the introduction of more specific boundary conditions precludes the resolution of the 
model in closed analytic form. 

On the other hand, when resorting to a numerical simulation of the observables one may take advantage 
of finite size effects to infer the critical behavior in the thermodynamic limit [Q. It happens, though, 
that finite size effects depend in general on the boundary conditions, in a way that may not be crucial 
but one cannot predict. There are again a few number of situations in which the asymptotic dependence 
on the spatial dimensions of the lattice has been rigorously studied. One of these cases corresponds to 
the analysis by Ferdinand and Fisher of the two-dimensional Ising model on large toroidal lattices 
The conclusions reached there support, in essence, the assumptions made in discussing finite size effects 
and, more precisely, the hypothesis of finite size scaling ^. Some open questions are raised, however, 
regarding the approach to the critical coupling, which is drastically influenced by the shape of the torus. 

The purpose of the present paper is to investigate finite size effects and critical properties of the Ising 
model on a class of two-dimensional lattices with spherical topology. Our choice for the elements of this 
class is not arbitrary, but is rather dictated by a prescription which makes possible increasing the size of 
the lattice without changing the local geometry. We propose to consider, in fact, a type of honeycomb 
lattices folded on the tetrahedron, which are built by assembling triangular blocks of the kind shown 
in Figure 1 as the faces of the polyhedron. One may construct a whole family of these lattices with 
increasing size, in such a way that the member of the 7V'^ generation A„ has a number of lattice points 
equal to n = 12N^. The coordination number is constant in each lattice. Moreover, from the point of 
view of the simplicial geometry, the curvature is always concentrated at the same faces, which are those 
formed by the three-fold rings around the four vertices. In principle, this makes the problem of taking 
the thermodynamic limit along our sequence of lattices well-defined. In rcf. Q] clear evidence was given 
of critical behavior in the ferromagnetic regime, as well as evidence supporting the hypothesis of finite 
size scaling applied to the model on the curved surface. 

The lattices we are considering may be understood as the result of applying nontrivial boundary 
conditions for the honeycomb lattice on the plane, though they have the effect of introducing curvature 
in the model. We undertake the investigation of the influence of these boundary conditions on flnite size 
effects and, more significantly from the physical point of view, on the critical properties of the model. 
Regarding the first point we will see that a discrepancy arises between the scaling of pseudocritical 
coupling constants for finite lattices and the true scaling of the correlation length. The second issue may 
be probably addressed in the continuum. In fact, investigations of the effect of boundary conditions in 
conformal field theories have been carried out before Anyhow, the inclusion of a conical singularity 
requires the kind of boundary condition which may call for a nonlocal operator in the theory |^ , so the 
analysis of our model in the continuum does not appear quite straightforward. 

The content of the paper is distributed as follows. In Section 2 we review the dimer approach applied 
to the computation of the partition functions and correlation functions of the model. Section 3 is devoted 
to the finite size analysis of data obtained with the above method, clarifying the issue concerning the ly 
critical exponent. In Section 4 we give technical details of the Monte Carlo simulations carried out to 
measure in some of the larger lattices. Section 5 contains the results for the critical exponents a, /?, 7 
obtained after combining data from the dimer approach and the Monte Carlo simulations. Finally we 
draw our conclusions and outline further directions of our work. 
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2 The dimer approach to the Ising model 



We review in this section the dimer formulation of the two-dimensional Ising model pO[ |T^. This 
approach presents the advantage of being applicable to lattices with an irregular coordination. It allows 
to write partition functions and correlation functions in closed compact form, in terms essentially of 
determinants of some coordination matrices for the lattice. This is something that cannot be achieved for 
our curved lattices by any other standard resolution method of the two-dimensional Ising model. There is 
no obvious way, for instance, as to how to apply the transfer matrix method to write down the partition 
function for a lattice with the topology of the sphere, not even to produce a numerical computation of 
the same. Within the dimer approach one may, in principle, compute the partition function for any of 
the hexagonal lattices inscribed on the tetrahedron. Although we have not been able to infer from this 
construction the thermodynamic limit along the sequence of growing lattices, the method shows very 
efficient to calculate observables like the specific heat or the correlation length with arbitrary precision. 
One can easily progress to lattices with up to more than 1000 points, with the possibility of applying a 
finite size analysis to study the critical behavior of the model. 

The dimer formulation of the Ising model makes use first of the equivalence between the partition 
function of the model and the dimer partition function of certain decorated lattice built from the original 
one One may apply afterward powerful techniques developed to perform the sum over dimer con- 
figurations. Let us review the former correspondence for the hexagonal lattice, while allowing for some 
kind of frustration which keeps constant the coordination number over the lattice. Given the collection 
of spins {cTi}, with i running over all the lattice points, the partition function Z is defined as the sum 
over all possible configurations 

where the Hamiltonian H is given by the sum over all the lattice links 

<jj> 

The factor \/kT is absorbed for simplicity in the definition of the nearest-neighbor coupling J. It is well 
known that ([^) can be cast as a sum over all the closed loops over the lattice. Calling this collection 
and {ni\ the respective numbers of links of the paths, we have actually 

Z = (cosh jy ^(tanh J)"' (3) 

/ being the total number of links of the lattice. One can draw a correspondence between each closed 
path and a dimer configuration in the appropriate decorated lattice. This is formed in our model by 
inserting a triangle in place of each of the points of the original lattice. To each of these one may assign 
four different states, depending on whether the point is traversed or not by a closed path and on what 
direction, in the first instance. These states are labeled in Figure 2. In a similar fashion, there are four 
different possible configurations of dimers on each triangle and adjacent links of the decorated lattice, 
which bear a one-to-one correspondence with the above states. These dimer configurations are labeled in 
Figure 3. One may easily convince that, establishing as a rule the equivalence between respective states 
in Figures 2 and 3, a unique closed path can be reconstructed starting from a given dimer configuration 
in the decorated lattice, and viceversa. Furthermore, if a weight equal to z = tanh J is given to each 
dimer on a triangle link and equal to 1 for dimers joining neighboring triangles, it is clear that the dimer 
partition function reproduces the sum in the Ising partition function (^). 

There exist, in turn, powerful techniques developed to the computation of dimer partition functions, 
which rely mainly on the relation between these and the Pfaffians of appropriate coordination matrices 
on the decorated lattice. We sketch here this relation, which has been worked out quite rigorously in ref. 
[ pl] |. The first step is to establish an order relation among the points of the decorated lattice. Once this 
is done, one may assign a matrix element Op^pj for points numbered pi and p2 such that ap^p^ — tanh J 
if the points belong to the same triangle, ap^p^ = 1 if they are nearest neighbors belonging to different 
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triangles, and ap-^p^ = otherwise. The sum over ah dimer configurations weighted as proposed before 
amounts to perform the sum over aU permutations {pi,p2, • ■ -Pk} Q 



'PlP2"P3P4 



(4) 



restricted by pi < < . . . pk-i and pi < p2, P3 < P4, ■ ■ ■ Pk-i < Pk- While there is no known algorithm 
to compute efficiently a sum of the above kind, one may think of allowing for an antisymmetric matrix 
A — {ttij}, so that the dimer partition function (a sum of all positive terms) may become proportional to 

with pi < p3 < . . pi < P2,P3 < P4, ■ ■ -Pk-i < Pk, as before, and (—1)^ being the signature of 

the permutation. The expression (H) reproduces the definition of the Pfaffian of the matrix A, which 
may be subsequently computed as the square root of its determinant. The remarkable conclusion which 
follows from the work of ref. |Q is that, for planar lattices, it is always possible to choose the sign of the 
nearest-neighbor matrix elements ap-^p^ so that all the terms in the sum (||) have the same sign. Since 
the matrix A becomes antisymmetric, it is customary to fix pictorically the sign of each ap-^p^ by giving 
an orientation to every link of the decorated lattice — cipipa is positive, for instance, if the arrow goes 
from pi to p2 — . We may enunciate the Kasteleyn theorem by saying that in any planar lattice there is 
always a system of arrows such that the dimer partition function can be computed as the Pfaffian of the 
corresponding antisymmetric coordination matrix. 

The lattices we consider here fall into the category of planar lattices since they have the topology of 
the sphere. As long as we are interested on dimers mainly for computational purposes, we simply give 
the recipes which have to be followed to form the appropriate system of arrows on a planar lattice. Once 
superposed on the plane, the lattice is made of so-called elementary polygons, which are the closed cycles 
that do not contain points in their interior. On the other hand, a polygon is said to be clockwise odd if 
the number of arrows pointing in the clockwise direction in the polygon is odd. The basic results which 
hold for planar lattices are that (a) it is always possible to choose a system of arrows such that all the 
elementary polygons are clockwise odd, and (b) with this choice and taking Up-^p^ as positive when the 
arrow goes from pi to p2, all the terms in the expansion of the Pfaffian have the same sign. 

In our case, a possible system of arrows realizing the property (a) for a decorated lattice inscribed 
on the tetrahedron is shown in Figure 4, where all the arrows for the triangle links (not drawn) are 
supposed to point in the clockwise direction. The advantage of this choice of arrows is that it keeps a 
regular pattern in the bulk, while only a few arrows on the boundary links have to be flipped to make all 
elementary polygons clockwise odd. In general, progressing to the next member of our family of lattices 
on the tetrahedron just amounts to adding a column of (decorated) hexagons at each side of Figure 4, 
expanding appropriately the vertical dimension. The system of arrows proposed may be extended in a 
straightforward way to larger lattices. According to the above discussion, we set always the absolute 
value of matrix elements Qp-^p^ equal to tanh J for points in the same triangle, and equal to 1 for nearest 
neighbors on different triangles. The partition function for any honeycomb lattice on the tetrahedron can 
be represented, therefore, in terms of the respective matrix A by 

Z= (coshJ)'(detA)i/2 (6) 

We have made use of the representation (H) to perform the numerical computation of the maximum 
of the speciflc heat (in the ferromagnetic regime) for lattices up to A1452. We have been able to measure 
that quantity with a relative error of less than 10"^ in most of the cases. Correspondingly, a precise 
determination of the coupling constant at which the maximum is attained in each lattice has been also 
possible (see Table 1). The values of these pseudocritical coupling constants are fundamental ingredients 
for the finite size analysis to be accomplished in the next section. We have also computed the values of 
the specific heat of the curved lattices at the critical coupling constant of the planar honeycomb lattice 
(see Table 2). These are also relevant under the hypothesis of finite size scaling since, as we will see, the 
sequence of pseudocritical temperatures converges in the thermodynamic limit to the critical temperature 
of the planar hexagonal lattice. 



is the total number of points in the decorated lattice. 
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We conclude this section with an outUne of how the two-point correlation functions can be obtained 
within the dimer approach ]lO| , [TH . Given two arbitrary spins dp and cTq in the lattice, the average 

WpUq) = 4 XI cTpCTge--^^ (7) 

may be computed with the following trick. One chooses a path C from Up to aq on the lattice, which will 
comprise a number of consecutive spins {(7p^ , (Tp^ , • ■ • }• The two-point function may also be expressed 
as ^ 

{<7pUq) — — CpCTpiCTpi . . . f7p„_lCTp,,^(Tqe~'^^ (8) 

^ <T.=±1 

Now we have that (JpCTp^, Op^Op^, . . .ap^^ciq are pairs of nearest-neighbor spins. Therefore, we find 

{(Tpaq) = — X Y[ (cosh J -I- (JiCTj sinh J) Y[ (sinh J -I- (7^ cr; cosh J) (9) 

where the first product extends to all the links which do not belong to C, and the second product runs 
over the links that do belong to the path. From (H) we arrive finally at the expression 

{Tpiq) = |(coshJ)'(tanh J)"+i|](tanhJ)"-'-'^^-l^ (10) 

where the sum, as in expression is over all the closed loops on the lattice, but with the difference now 
that the number of links in each loop belonging to C have to be weighted with (tanh J)^^ rather than 
with tanh J. It becomes obvious that all the machinery of the dimer formulation can be applied again to 
transform the sum in ([l^) into a suitable dimer partition function on the decorated lattice, so that 

(apaq) = l(cosh J)'(tanh J)'"+i(det A')^/^ (11) 

The appropriate coordination matrix A' — {a'^-j has to keep track of the different weight that the rn + 1 
links in C carry in the sum over closed loops. 



3 Finite size scaling and critical exponents 
3.1 Finite size scaling 

It is well known that singularities in the free energy (i.e. phase transitions) can only occur in the 
thermodynamic limit. For finite volumes the free energy is an analytic function of the temperature and 
any other parameter in the Hamiltonian. The thermodynamic singularities are thus smoothed out around 
the transition point. A trace of the existence of such non-analyticities is the presence of some peaks in the 
specific heat CV or magnetic susceptibility x curves. The dependence on the linear size of the system L of 
the location of the maxima of those peaks and their height permits the description of the thermodynamic 
limit from finite-size data [Q. 

In second order phase transitions this round-off is due to the fact that the correlation length ^ is 
limited by the size of the system. This fact defines a pseudocritical coupling J*{L) such that 

ar{L)) ^ L (12) 

At this point the surface contribution to the free energy is not negligible compared to the bulk one. In 
the vicinity of a second order transition point Jc the correlation length diverges with a power-law given 

by 

^J) - {J - Jc)-" (13) 
From ( p^ and (|l^) it can be derived the dependence of the pseudocritical coupling on the lattice size 

\J*iL) - J,| ^ L-i/" (14) 
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Unfortunately in practical situations it is a very involved task to compute such quantity J*{L). It is 
easier to look at the position and height of those peaks mentioned above. If some quantity P diverges 
near the critical point as ^ 

P{J)^\J-Jc\-''; p>0 (15) 

then it can be shown Q that for a finite volume it attains a maximum value Pmax{.JL) at a point JL{Pmax) 
given by 

UdP^ax) - Jc\ - L-^" (16a) 

Pmax{JL) - L"!- (16b) 



when L is large enough. In most systems it is found that 

V 



(17) 



but this is not a general result. There are some examples where this property does not hold: the spherical 
model, the ideal Bose gas p2| and the one-dimensional q = oo clock model |T3| . In the present paper we 



will face another situation in which the relation ( |17|) is violated. On the other hand, if 6* > l/i^ then the 



behavior of this quantity at finite volume evaluated at the critical point Pl{Jc) is the same as in (16b) 

Pl{Jc)-Lp'- (18) 



Using (16a) and (16b) the critical coupling Jc and the critical exponents ratio p/v can be derived from 



finite-size data. When Jc is explicitly known and Op > 1/v then (18) can be used instead of (16b). In 
this paper we are mainly concerned with the analysis of the susceptibility and the specific heat. So, we 
will obtain estimates of Jc and the ratios 7/1/ and a/i^. The rest of the critical exponents may be derived 
using the scaling relations |Q : 

(19a) 
(19b) 
(19c) 
(19d) 





1 
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V 
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4 


- 1 



However, in this paper we will check numerically scaling relations (19a) and (19b). We will obtain 
independent estimates of v and fijv respectively in terms of the analysis of the correlation length (see 
below) and the magnetization at the critical point (see Section 5). 



3.2 Critical point and exponent Qcv 

In Section 2 we showed that for lattices up to A1452 we were able to obtain very accurate estimates of 
the internal energy E and the specific heat. These quantities are defined hereafter as follows 

E« = ^iH) (20) 

where the factor 3F/2 is equal to the number of links in a lattice of V sites and (t(-) is the standard 
deviation. Thus, the values of CVmax and JhiCVmax) can be computed with high precision. In what 

^Hereafter we will denote quantities computed in a finite volume with a subscript L meaning the linear 
size of the system (i.e. Pl)- Whenever no subscript is present, the thermodynamic limit is assumed (i.e. 

P = limL^ooPL). 
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follows we will identify the lattice linear size L with the index N characterizing the fuUerene lattice (see 
Section 1). This choice is consistent as the volume increases as = 12 N"^ . 
Data will be fitted according to Equation (16a) 

JN^CV^ax) ^Jc + AN-'^''^ (22) 

using a least method. Here the input errors are given by the precision of the computer in calculating 
JNiCVmax)- III order to obtain a more reliable result, we will sequentially remove the point with smallest 
A''. One eventually can observe a monotonous trend to some value, which will be identified with the 
thermodynamic limit. 
Our best result is 

Jc = 0.65850 ±0.00002 (23) 

for 5 < iV < 11. In this case = 0-07 with 2 degrees of freedom. Throughout this paper all the errors 
associated with our final results will be equal to 2 standard deviations — i.e. 95% of confidence level — . 
The later result is compatible with the critical point of the Ising model on a toroidal honeycomb lattice 

Jtorus ^ 1 ^Qg(2 + V3) = 0.65848 (24) 

This result can be easily derived using duality ||lj, |l5|. Thus, our data strongly supports that both 
critical points coincide. 

A good estimate of 9cv is obtained by repeating the fit with Jc = J'°'"^. The result is 

6icy = 1.745 ±0.015 (25) 

for 6 < A' < 11 and with = 1-2. This one is in clear disagreement with the result expected for a lattice 
on a torus {9cv = 1 (^)- This fact makes necessary a direct determination of the v critical exponent, in 
order to see if the above measurement bears any relation to it (see below). 

We will skip here the analysis of CVmax, as it is rather subtle to distinguish between a logarithmic 
and a power-law behaviors, specially when such a power is rather small. This will carried out in Section 
5. 



3.3 Correlation length and exponent ly 

An independent way to compute the critical exponent ly is to study the correlation length near the critical 
point. This quantity is defined in terms of the connected two-point function 

(o-Q ■ (^r)" = (c^o ■ O-r) - (cto) (Cr) ^ 6^''/^ (26) 

when r is large enough. The connected two-point correlators are equal to the usual ones (uo • CTf) for 
J < Jc (unbroken phase). This feature allows their exact computation using the machinery developed in 
Section 2. (For finite lattices odd quantities such as the magnetization are always equal to zero, even in 
the broken phase). 

A major problem is how to recover ^Ar(J) from the finite- volume correlators {(Jq ■ o'r)jv For a torus 
of linear size L these functions are expected to behave as ~ cosh((a; — L/2)/S_l) when a; ^ 1. But this 
question is not clear for the truncated tetrahedron. On the other hand, it is well known [ pl| that the 
correlation length does depend on the direction along the spins Ur are disposed. However, the same 
critical behavior is expected for all the possible directions. 

This study has been carried out on the lattice Ag72 , which is the largest one allowed by our computer 
facilities. We believe that this one is large enough to see the thermodynamic limit. We have chosen 
couples of spins along the diagonal in the representation on Figure 4 of the tetrahedron unfolded on the 
plane. Our choice for a-r allowed us to introduce an increasing distance between spins running from 3,6,- • • 
up to 27 (in units of the lattice spacing). At r = 27 both spins are located at antipodal points. If we 
go on along the diagonal we finally arrive at ctq. For that reason we expect that the correlators behave 
for large r as a symmetrized version of equation (|2^). In order to improve our results, we include in our 
ansatz the correct leading term for the square lattice on the torus ||ll| 

(ao(Tr)-^e-^/«(l + 0(l/r)) (27) 
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suitably symmetrized around r = 27. 

We have analyzed the cases J = 0.58, 0.59, 0.60, 0.61 and 0.62 (see Table 3 and Figure 5). We have 
obtained extremely good fits for all these cases, giving differences of order ^ 10"^. Although the lattice 
A972, as any of the lattices inscribed on the tetrahedron, is not homogeneous, it is remarkable that the 
values of the two-point functions at each different J fit, to a high degree of precision, to the correct 
leading behavior for the Ising model on a square lattice on the torus. The deviation that we have found 
from the dependence (j2^) appears to be even smaller than for similar measurements carried out for the 
lattice on a torus. The estimated values of the correlation length are shown in Table 3 and will be used in 
the computation of the v critical exponent. We unsuccessfully tried to fit data for J > 0.62 to (|27|). The 
reason is that very close to the critical point we have to take into account in Equation ( p7| ) the 0(l/r) 
(or even higher) terms. 

Given the values of ^jvLgl"^) Table 3, we tried to fit them according to ( p^ ) 

CN=<,{J)=A\3-3.r (28) 

We obtain a value equal to v = 1.01 ± 0.04 with = 0.36. However, if we drop the point with J = 0.62 
(the closest to Jc) the result is 

= 1.00 ±0.06 (29) 

with a remarkable small value of the ^ 7 • 10^^. 

The conclusions of the analysis of the data coming from the dimer computations can be summarized 
in the following points (a) The critical point is compatible with J*"''^**. (b) The critical exponent v — 
1.00 ± 0.06. (c) The finite-size exponent dcv — 1.745 ± 0.015 is significantly different from In 
summary, our results suggest that the critical properties of the Ising model on the truncated tetrahedron 
are the same as on the torus. However, we find a very clear difference in the finite-size behavior of those 
models, as long as the scaling of pseudocritical coupling constants (determined from the maxima of the 
specific heat) does not match with the scaling behavior of the correlation length. On the tetrahedron the 
thermodynamic limit is achieved much faster than on the torus, at least in what concerns to the specific 
heat. 



4 Technical aspects of the Monte Carlo simulations 

We have performed several Monte Carlo (MC) runs for different lattice volumes V — 12N'^ and coupling 
constants J. The relevant information about the simulations can be found in Table 4. 

We have used a Metropolis algorithm with the R250 pseudorandom-number generator -initialized 
with the RANMAR subroutine-. The period of such generator is equal to 2^^° — 1. Recently, it has been 
claimed 
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that the combination of the Metropolis algorithm with the R250 generator gives better 
results than other more involved procedures. We have compared those values obtained both by the dimer 
approach and by direct MC simulation. They were consistent within statistical errors. 



In all cases we have measured the internal energy density and the magnetization defined as 



M 



N 



(30) 



We have also measured the specific heat and the magnetic pure-phase susceptibility, 



Xn = Va 



(31) 



In all cases we discarded the first 10^ MC steps for termalization. Then we have measured each 
observable once every typically 100 MC steps. In this way we obtained statistically independent data, as 
it can be checked by computing the corresponding autocorrelation times jl^ . 

In order to calculate de maximum value of the specific heat and the susceptibility we have used the 
Spectral Density Method ||l8|. At a given coupling J we can obtain the histograms N{E,M; J) which 
keep track of the numbers of configurations with magnetization M and internal density energy E. This 
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information is enough in order to compute the expectation value of any function f[M,E) at any other 
coupUng J'. In our case the magnetic field is zero and in the equations ( pO| , pTpO|j3l| ) the observables do 
not depend on E and M simultaneously. For those reasons we could use the following formulae 

lf(^,^^lr^ EEfimo{E;J)eM{J'-J)E} 

{J[E))[J ) = , , T--r rr— [62) 

T,E^o{E; J) exp{(J' - J)E} 

2^E^o(E;J)exp{( J' ~ J)E\ 
where the one dimensional histograms are defined as follows 

No{E;J) - ^iV(i;,M;J) (34) 

Nf{E-J) - ^/(M)iV(i?,Af; J) (35) 

M 

The Spectral Density Method gives the correct answer for couplings J' close to the coupling J where 
the simulation was performed. A criterion for the applicability of such method is the following 

In most cases only one simulation at J — Jc is enough in order to determine the maximum of CV and 
X- However, for the smaller lattices an additional run had to be performed in order to obtain a reliable 
estimate of such quantities. 

We have divided up the entire sample into typically 30-120 subsamples, each of them containing 
1000 measures. For each subsample we computed every quantity (including CVmax, JN{CVmax)i ■ ■ •)■ 
With these estimates we calculated the statistical errors using the jack-knife method ||l^. In this way, 
the effect of correlation among data was taken into account. 

We have done all the MC simulations on a VAX 9000 machine with a vectorial processor. The program 
was not fully vectorizable, as the lattice could not be splitted into two disjoint sublattices, in such a way 
every element of one sublattice is surrounded by elements belonging to the other one|^ However, we could 
divide the whole lattice into three subsets. The elements of the first two are arranged on two disjoint 
triangular sublattices, so their update could be fully vectorize. On the other hand, the rest of the spins 
can be located on two lines joining pairs of vertices on the tetrahedron and their number depend explicitly 
on the planar representation of the lattice. For these ones the update is clearly not vectorizable. However, 
their effect on the CPU time is not very important for the larger lattices, as their number behaves as 



5 Results of the Monte Carlo simulations 
5.1 Position of the critical point 

Here we will repeat the analysis of Section 3, but with all the data of Table 1. For data coming from the 
dimcr analysis the input error will be taken as the precision of the subroutines used. For those coming 
from the MC simulations the error is given by the jack-knife method described in the preceding section. 
Data for iV > 11 posses large error bars compared with the rest. For that reason the fits presented in 
Section 3 cannot have large variations. Our final results for the specific heat are 

Jc = 0.65850 ±0.00002 (37) 

for 5 < AT < 21 and with = 0.67. And 

61^/ = 1.745 ±0.015 (38) 
■^This feature has to do with the onset of frustration in the antiferromagnetic regime. 
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for 6 < TV < 21 and with = 1-9 (See Figure 7). 

If we repeat the same procedure with the susceptibihty we obtain an estimate for Jc compatible with 
the later one, but with larger error bar. If we fix this quantity to J^:™^^ we arrive at the following estimate 
for the exponent 9^. 

6'x = 1.01 ± 0.02 (39) 

for 9 < TV < 21 and = 0.8 (See Figure 7). 

We observe that 9^ is close the value 1/;^ = 1, in agreement with the Ising model on the torus. We 
conclude that 

\ = h< ^cv (40) 

and this feature implies that one can obtain the critical exponents using either eq. (16b) or eq. (|l8|). As 
we have identified the critical coupling of our system, it seems more natural to rely our conclusions on 
eq. (|l8|). In any case, the values obtained from eq. (16b) are always consistent with those presented in 
this paper within statistical errors. 



5.2 Exponent ratios j/u and /S/u 



The value of j/i^ can be derived using the value of the critical magnetic susceptibility xn{Jc)- We have 
fitted our data to 

Xn{Jc) = AN^/-- (41) 

and our best result is 

^ = 1.73 ± 0.02 (42) 

for values of N ranging from 21 to 9 and with x^ = 0.79 This result is very close to the usual Ising model 
one j/i^ = 7/4 = 1.75. If we fix j/i' to this value we obtain a x^ value of ~ 2.8 which shows that the fit 

is reasonably good (See Figure 8). 

To determine the value of f]/^ in an independent way (i.e. not using scaling relation (19a)) we will fix 
our attention on the magnetization at the critical point. In the thermodynamic limit the magnetization 
near Jc behaves as 

M{J)^\J-Jcf (43) 
Using similar arguments as in Section 3 it can be predicted that 

Mn{Jc) ~ N-"/'' (44) 

The result of performing such fit is the following 

- =0.126 ±0.004 (45) 
v 

using data with 7 < iV < 21 and with x^ = 0.33 (See Figure 9). This value is also very close to the Ising 
one P/v — 1/8 = 0.125. This result supports that the scaling relation ( |l9a| ) does hold in this model. 

This feature can be used in order to obtain a more accurate estimate of those exponents. We can try 
to fit xn{Jc) and M^^Jc) simultaneously using explicitly the relation ( |19a| ). The result is 

7/i/= 1.748 ±0.008 

/3/z^ = 0.126 ±0.004 ^ ' 

Thus, our results strongly suggest that the ratio j/v = 7/4 as in the Isi ng m odel on a torus. Notice 
that the error in that ratio is less than 0.6%. Using the relations (19c) and (19d) we can derive the value 
of the exponents ry and 5. 

r] = 2 - - = 0.252 ± 0.008 (47a) 

V 

6 = --1 = 14.8 ±0.6 (47b) 
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5.3 Exponent ratio a/u 

We can perform the same game for the specific heat and try to obtain the exponent ratio a /v. If we try 
to fit data to the function 

CVn[Jc) = A + BN""/" (48) 

we do not obtain a satisfactory result. The best one gives a ratio a/v ^ 0.060 with ^ 9 and 9 < iV < 21 
(See Figure 8). On the other hand, and motivated in part for the preceding results we could try to fit 
the data to a logarithmic function. 

CVn{Jc) = A + BXogN (49) 

In this case, the fit is successful giving a = 1-9 with 7 < < 21 (See Figure 10). This immediately 
implies that 

a = (50) 

and using equation (19b) 

v^l (51) 

which is in agreement with the result from direct measurements of the correlation length displayed in 
Section 3. Thus, both exponents take the same values as in the Ising model on the torus. On the other 
hand, we have also verified that the scaling relation (19b) does hold in this model. 



6 Conclusions and outlook 

In this paper we have presented the first deep study of the critical properties of an Ising model on a 
lattice with the topology of the sphere. In particular we have chosen the family of honeycomb lattices 
that can be constructed on the tetrahedron. Our results can be summarized as follows 



The dimer approach is very useful and competitive for lattices up to ^ 10'^ points. It provides very 
accurate data for the internal energy, specific heat and two-point correlators. 

The critical properties of the Ising model on the tetrahedron are just the same as on the torus. In 
particular we have checked that Jc is the same, as well as the critical exponents a, 7 and /3. 



• We also have checked two scaling relations (19a) and (19b) among these exponents. Using the two 
other equations we obtained the last two critical exponents 5 and 77. They are in agreement with 
those corresponding to the Ising model on a torus. 

• But the finite-size scaling properties of those two systems are not the same. In our case the position 
of the maxima of the specific heat scales near Jc with a critical exponent that does not bear any 
relation to the critical behavior of the correlation length. That is, the thermodynamic limit is 
achieved faster on the tetrahedron. However, the behavior of the susceptibility is just the same in 
both ones. 

Our results suggest that the same analysis carried out in other types of lattices with the same topology 
will yield the same conclusions: the critical behavior does not change, although variations in their finite- 
size properties are expected to hold. All of them belong to the same universality class of the Ising model 
on a torus. 

On the other hand, the Ising model on the tetrahedron can be view as the Ising model with some 
non-standard boundary conditions. These ones has the advantage that the critical behavior is reached 
before than for periodic boundary conditions. 

Ferdinand and Fisher |j studied how the exponent 9cv varies for the Ising model defined on a 
square lattice on a m x n torus. They concluded that Ji{CVmax) behaves as 

:Mg^ = i + M + o(n (52) 
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where / = (m^^ + n^^)^^ measures the hnear size of the torus and rj — m/n its shape. They also 
showed that 6(77) is not a monotonous function of 77. In particular 6(77) > in the range 77 G {Vo^tVo) 
with 770 — 3.139278, and 6(77) < in 77 G (0,7;q"^) U (770, cxd). Exactly at 7/ = 770 and 77q"^ the function 
6(7/) vanishe^, so at these points the leading term in (52) vanishes and its behavior is controlled by the 
subleading term. In this model it can be written as 



Jl{CVmax ) 



clog I 

1 - + ■■■ (53) 



The behavior of 6(77) as a function of the shape of the torus is explained as a highly non-trivial interplay 
between the different terms which appear in the expression of the partition function. We believe that 
the same feature is present in the Ising model on the truncated tetrahedron. In this case, the shape of 
the lattice is fixed, but the chosen boundary conditions are the basic ingredient which makes the leading 
term in ( js^ ) vanish. 

We should mention that the our conclusions may not apply to the antiferromagnetic regime. The 
reason is that for such boundary conditions, the lattice is not bipartite. Thus, the phenomenon of 
frustration may occur in that regime. This feature is absent in the Ising model on the torus. In this case, 
the lattice is bipartite and for sufficiently low temperatures we find a Neel ground state. This question 
on the tetrahedron is currently under research. 

Finally, we would like to say a few words about the continuum theory which is attained in the 
thermodynamic limit. The evidence we have found of scaling suggests a description in terms of the 
fields and weights of a conformal field theory. This cannot be a trivial example of field theory on the 
sphere (or on the plane), since four curvature singularities arise which cannot be removed by conformal 
transformations. Under the assumption of conformal invariance, though, we could still stick three of the 
vertices together (at a point we may take as infinity) by means of SL(2, C) transformations, leaving alone 
a singularity in the bulk. This picture is close in spirit to the Coulomb gas representation of conformal 
field theories, but with the difference now that not all the curvature is pinched at the point at infinity. A 
conical singularity on the plane may have sensible effects on the correlation functions of the theory. We 
recall here another example with nontrivial boundary conditions, namely that of a conformal field theory 
on the semiplane. In this case correlators which are taken at a finite distance of the boundary do not 
measure the conformal weights of the theory on the plane pO[ | . In our model, it may not be necessary to 
compute correlators infinitely far away from the singularity to measure bulk conformal weights, but again 
some dependence on the location of the points should be expected. This point should deserve further 
clarification, though its investigation in the lattice would require more powerful computer facilities than 
used in the present work. 
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Figure Captions 

Figure 1: Generic triangular block for honeycomb lattices embeded on the tetrahedron. 

Figure 2: Different paths through a honeycomb lattice site. 

Figure 3: Dimer configurations around a triangle of the decorated lattice. 

Figure 4: Decorated lattice for the second generation. The outer lines show the identifications of 
boundary links which embed the lattice on the tetrahedron. 

Figure 5: Values of the unsubstracted two point correlators (cTrCJo) for different separations r and 
couplings constants J = 0.57, 0.58, • • •, 0.62. The least fits are also despicted. 

Figure 6: Values of ^{J)~^ for those J listed in Table 3. The straight line corresponds to the fit. 

Figure 7: Values of the position of the maxima of the specific heat (circles) and the magnetic suscepti- 
bility (squares) with respect to the critical coupling J*°''"^. The power-law fits are also despicted. 

Figure 8: Power-law fits of the values of the specific heat (squares) and the magnetic susceptibility 
(circles) at the critical point Jc- 

Figure 9: The same as in Figure 8 for the magnetization Mm{Jc)- 
Figure 10: Logarithmic least fit to CVn{Jc)- 
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Table Captions 



Table 1: For each lattice size N we show the values of the maximum of the specific heat and the position 
of such point {CVmax, JN{CVmax)) and the same for the susceptibility {xmax, JNiXmax))- Those values 
marked with a f were computed using the exact partition function 0|, those with a % using the dimer 
approach (see Section 2) and those with an ★ by evaluating exactly (and numerically) the partition 
function. The rest were obtained by means of the Monte Carlo simulations described in Section 3. 

Table 2: For each lattice size N we show the values of the specific heat CVn{Jc)i the susceptibility 
Xn{Jc) and the magnetization Mn{Jc) evaluated at the critical point J*°''"''. The symbols posses the 
same meaning as in Table 1. 

Table 3: Values of the inverse of the correlation length computed for the = 9 lattice and different 
values of the coupling J. 

Table 4: Number of Monte Carlo steps performed for each simulation in a lattice of volume V = 12N^ 
and coupling constant J. 
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TV 




JN{CVmax) 


X,max 


JNiXmax) 


1 


1.43923551(l)t 


0.467332(l)t 


1.2564278(1)* 


0.487805(1)* 


2 


1.65204595(l)t 


0.608224(1)^ 






3 


1.89238232(l)t 


0.634238(l)t 


7.91(2) 


0.6215(3) 


4 


2.07977738(1)* 


0.643862(1)* 






5 


2.23081551(1)* 


0.648580(1)* 


18.98(7) 


0.6377(3) 


6 


2.35682745(1)* 


0.651275(1)* 






7 


2.46477396(1)* 


0.652973(1)* 


34.4(1) 


0.6444(2) 


9 


2.64288048(1)* 


0.654929(1)* 


53.3(2) 


0.6480(2) 


11 


2.786536(6)* 


0.65598(5)* 


75.9(4) 


0.6498(2) 


15 


3.02(1) 


0.6571(3) 


128.4(8) 


0.6523(2) 


21 


3.30(2) 


0.6579(3) 


232(2) 


0.6540(2) 



Table 1 



N 


CVn{Jc) 


XJv(Jc) 


Mn{Jc) 


1 


0.99350489(1)* 


0.8617020(1)* 


0.82868806(1)* 


2 


1.51639951(1)* 






3 


1.81746516(1)* 


6.889(3) 


0.6904(9) 


5 


2.19465654(1)* 


16.89(9) 


0.645(1) 


6 


2.32898639(1)* 






7 


2.44249002(1)* 


30.8(1) 


0.6168(7) 


9 


2.62744012(1)* 


48.2(3) 


0.597(1) 


11 


2.77(1) 


68.3(4) 


0.5828(9) 


15 


2.998(9) 


116(1) 


0.560(1) 


21 


3.24(1) 


210(2) 


0.536(2) 



Table 2 
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J 


^JV=9 


0.57 


0.204(7) 


0.58 


0.181(5) 


0.59 


0.158(1) 


0.60 


0.135(2) 


0.61 


0.112(1) 


0.62 


0.088(1) 



Table 3 



N 


J 


MC steps 


3 


Ttorus 


3 • 10*^ 


3 


0.62 


3 • 10*^' 


5 


^torus 


3 • 10^ 


5 


0.64 


3 • 10*^ 


7 


Ttorus 
•^c 


12 • 10*^' 


9 


Ttorus 
•^c 


12 • 10*^ 


11 


Ttorus 
•^c 


12 • 10*^ 


15 


Ttorus 
•^c 


11 • 10*^ 


21 


Ttorus 


12 • 10^ 



Table 4 
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Figure 1: 
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